package br.unb.statistic.distribution;

//Classes necessárias
import br.unb.statistic.distribution.BetaDistribution;
import br.unb.statistic.distribution.NormalDistribution;
import br.unb.utils.Resource;
import java.io.BufferedWriter;
import java.io.IOException;
import java.text.SimpleDateFormat;
import java.util.Date;

/**
 *
 * Gerador de tráfego multifractal utilizando o algoritmo MMW (Modelo Multifractal Baseado em Wavelets) criado por Flávio Henrique Teles Vieira
 * em sua tese de doutorado 'Contribuições ao Cálculo de Banda e de Probabilidade
 * para Tráfego Multifractal de Redes', Campinas São Paulo, 2006
 *
 * Data: 15 de agosto de 2007
 * @author Tiago Fraga
 *
 * Multifractal traffic generator using the MMW algorithm (Wavelet-Based Multifractal Model)
 * created by Flávio Henrique Teles Vieira
 * In his doctoral thesis "Contributions to the Calculus of Probability and Bandwidth
 * For Multifractal Network Traffic ", Campinas Sao Paulo, 2006

 * Data: August 15, 2007
 * @author Tiago Fraga
 */
public class MMWGenerator {

    // Declaration of input parameters of the algorithm MMW
    int n; // Number of samples of traffic to be generated
    double alpha; //alpha parameter
    double gamma; //gamma parameter
    double rho; //rho parameter

    /**
     * Constructor of MMW traffic generator, with default parameters
     * Default constructor uses size 4 and parameters set to 1
     */
    public MMWGenerator() {

        //"Passo 1 - Configuração Inicial: Conjunto de parâmetros de entrada (alpha, gamma, rho) para uma série de tráfego de tamanho 2 elevado a n amostras;"
        //"Step 1 - Initial Setup: Set of input parameters (alpha, gamma, rho) for a series of traffic size (2 power 2) = 4 samples;"
        this.n = (int) Math.pow(2, 2);
        this.alpha = 1D;
        this.gamma = 1D;
        this.rho = 1D;
    }

    /**
     * Constructor of MMW traffic generator, with custom parameters
     *
     * @param n - Number of samples of traffic to be generated
     * @param alpha
     * @param gamma
     * @param rho
     */
    public MMWGenerator(int n, double alpha, double gamma, double rho) {

        //"Passo 1 - Configuração Inicial: Conjunto de parâmetros de entrada (alpha, gamma, rho) para uma série de tráfego de tamanho 2 elevado a n amostras;"
        //"Step 1 - Initial Setup: Set of input parameters (alpha, gamma, rho) for a series of traffic size (2 power 2) = 4 samples;"
        this.n = n;
        this.alpha = alpha;
        this.gamma = gamma;
        this.rho = rho;
    }

    /**
     * Method that generates a lot of traffic with multifractal distribution
     * according to the synthesis algorithm of traffic MMW.
     *
     */
    public double[] getMMWTraffic() {

        //"Passo 2 - Inicie com a escala j = 1, correspondendo é fração de agregação igual a 2 elevado a j = 2;"
        //"Step 2 - Start with the scale j = 1, corresponding aggregate fraction is equal to 2 raised to j = 2;"
        int j = 1; //Variable-scale representation

        double[] processVariance = new double[n]; // Storage array for variances of scale
        double[] processZ = new double[n]; // Storage array for Z variable
        double[] waveletCoefficientSecMoment = new double[n]; // Storage array of second-order moment of wavelet coefficients
        double[] traffic; //Array containing the generated traffic

        do {
            //"Passo 3 - Calcule a variáncia do processo agregado var[X elevado a m] onde m = 2 elevado a por meio da equação '3.12' ( getProcessVariance() )"
            //"Step 3 - Calculate the variance of the aggregate process var [X high am] where m = 2 ^ from the equation '3 .12 '(getProcessVariance ())"
            processVariance[j] = this.getProcessVariance(j);

            //"Passo 4 - Calcule Zj usando a variância do processo agregado var[X elevado a m] calculado no passo anterior utilizando a equação '3.19' ( getProcessZ() )"
            //"Step 4 - Calculate ZJ using the variance of the aggregate process var [X high am] calculated in the previous step using the equation '3 .19 '(getProcessZ ())"
            processZ[j] = this.getProcessZ(processVariance[j], j);

            //"Passo 5 - Estime o momento de segunda ordem dos coeficientes de wavelet E(w2) do processo X(t) via equação '3.18' ( getSecOrderMomentConfWavelet() ) e pelo valor de Z calculado no passo anterior"
            //"Step 5 - Estimate the second order moment of wavelet coefficients E (w2) of the process X (t) via equation '3 .18 '(getSecOrderMomentConfWavelet ()) and the value of Z calculated in the previous step"
            waveletCoefficientSecMoment[j] = this.getConfWavSecOrderMoment(processZ[j]);

            //"Passo 6a - Incremente j de 1. Se j != n, vá para o passo 3"
            //"Step 6a - Increment j, 1. If j <> n, go to step 3"
            j = j + 1;

        } while (j != n);

        //"Passo 6b - Se j == n, calcule a média e a variância dos coeficientes em escala de escala mais grosseira (U0,0) usando as equações '3.20' e '3.21'"
        //"Step 6b - If j == n, calculate the mean and variance of coefficients on a scale from coarser scale (U0, 0) using equations '3 .20 'and '3 .21'
        processVariance[0] = this.getProcessVariance(0);
        processZ[0] = this.getProcessZ(processVariance[0], 0);
        //double meanCoarsetScaleCoefficient = Math.pow(2, n/2)*(Math.exp(rho+Math.pow(gamma, 2)/2));
        double meanCoarsetScaleCoefficient = (Math.exp(rho + Math.pow(gamma, 2) / 2));
        double varCoarsetScaleCoefficient = (Math.exp(4 * rho + 4 * Math.pow(gamma, 2)) * Math.pow((alpha + 1 / alpha + 0.5), 2)) + 2 * processZ[0] - (Math.exp(2 * rho + Math.pow(gamma, 2)));

        //"Passo 7 - Calcule os valores de pj usando a equação (3.34)"
        //"Step 7 - Calculate the values of pj using equation (3.34)"
        double[] p = new double[n];

        //Calculando p0 usando a fórmula (2.68)
        //Calculating p0 using the formula (2.68)
        waveletCoefficientSecMoment[0] = this.getConfWavSecOrderMoment(this.getProcessZ(this.getProcessVariance(0), 0));
        p[0] = ((varCoarsetScaleCoefficient / waveletCoefficientSecMoment[0]) - 1) / 2;

        j = 1;

        do {
            p[j] = (((waveletCoefficientSecMoment[j - 1] * (p[j - 1] + 1)) / waveletCoefficientSecMoment[j]) - 1) / 2;

            j = j + 1;

        } while (j != n);

        //"Passo 8 - Aqui começa um procedimento recursivo tendo obtido os valores de pj's. Mude j=0 e calcule os coeficientes de escala na escala mais grosseira (U0,0)"
        //"Step 8 - Here begins a recursive procedure and obtained the values of pj's. Change j = 0 and calculate the scale coefficients at coarser scale (U0, 0)"
        j = 0;
        double[] scaleCoefficient = new double[1];

        //Calculo de U0,0 utilizando a distribuição Normal;
        //Estimated U0, 0 using the Normal distribution
        NormalDistribution gaussian = new NormalDistribution(meanCoarsetScaleCoefficient, varCoarsetScaleCoefficient);
        scaleCoefficient[0] = gaussian.rand();

        //"Passo 9 - Gere amostras aleatórias correspondentes a Aj,k para k = 0, 1,...(2 elevado a j) - 1
        //"Step 9 - Generate random samples corresponding to Aj, k for k = 0, 1, 2 ,...( high j) - 1
        BetaDistribution beta;
        double[] multiplierA;
        do {
            beta = new BetaDistribution(p[j], p[j], -1, 1);
            multiplierA = new double[scaleCoefficient.length];

            for (int i = 0; i < multiplierA.length; i++) {
                multiplierA[i] = beta.rand();
            }

            scaleCoefficient = this.getNextScaleCoefficients(scaleCoefficient, multiplierA);

            //"Passo 11 - incremete j de 1. Se j < N, repita os passos 9 e 10
            //"Step 11 - incremete j 1. If j <N, repeat steps 9 and 10
            j = j + 1;

        } while (j < n);

        //"Passo 12 - A série de tráfego gerada é X(k) = (2 elevado a -j/2)*Uj,k , k = 0, 1,..., 2 elevado a N - 1"
        //"Step 12 - A series of traffic generated is X (k) = (2 ^-j / 2) * Uj, k, k = 0, 1 ,..., 2 ^ N - 1"
        traffic = new double[(int) Math.pow(2, n)];
        for (int i = 0; i < traffic.length; i++) {
            traffic[i] = Math.pow(2, -j / 2) * scaleCoefficient[i];
        }

        return traffic;
    }

    /**
     * Variance of an aggregate process X j
     *
     * @param j - Scale of agregated process
     * @return varP - variance of process
     */
    private double getProcessVariance(int j) {

        double varP;
        double part1, part2, part3;

        //Fórmula de cálculo
        part1 = Math.exp(2 * rho + 2 * Math.pow(gamma, 2));
        part2 = Math.pow((alpha + 1) / (alpha + 0.5), n - j);
        part3 = (Math.exp(2 * rho + Math.pow(gamma, 2))) * Math.pow(2, 2 * j - 2 * n);

        varP = part1 * part2 - part3;

        return varP;
    }

    /**
     * Retorns the Z variable of a given process with scale 'j' and variance 'var'
     *
     * @param var - Process Variance
     * @param j - Scale of agregated process
     * @return z - Z variable of the process
     */
    private double getProcessZ(double var, int j) {

        double z;
        double part1, part2, part3;

        part1 = Math.exp(2 * rho + 2 * (Math.pow(gamma, 2)));
        part2 = (alpha + 1) / (alpha + 0.5);
        part3 = Math.exp(2 * rho + (Math.pow(gamma, 2)));

        z = (Math.pow(2, j - 1)) * (var - part1 * part2 + part3);

        return z;
    }

    /**
     * Calculation of second order moments of wavelet coefficients
     *
     * @param z - Z variable
     * @return secMoment - Second order moment of the process
     */
    private double getConfWavSecOrderMoment(double z) {

        double secMoment;
        double part1, part2;

        part1 = Math.exp(4 * rho + 4 * (Math.pow(gamma, 2)));
        part2 = Math.pow((alpha + 1) / (alpha + 0.5), 2);

        secMoment = part1 * part2 - 2 * z;

        return secMoment;
    }

    /**
     * Calculates the coefficients of scale finer scale (close range)
     *
     * @param lastScaleCoefficients - coefficients from previous scales
     * @param mulA - Multipliers A
     * @return
     */
    private double[] getNextScaleCoefficients(double[] lastScaleCoefficients, double[] mulA) {

        int k = lastScaleCoefficients.length; //current scale
        double[] nextScaleCoefficents = new double[2 * k];
        double waveletCoefficient;

        for (int i = 0; i < k; i++) {

            //"Passo 9 - Calcule os coeficientes de wavelet na escala j (Wj,k) usando (2.66)"
            //"Step 9 - Calculate the wavelet coefficients at scale j (Wj, k) using (2.66)"
            waveletCoefficient = lastScaleCoefficients[i] * mulA[i];

            //"Passo 10 - Na escala j, calcule os coeficientes de escala na escala j + 1 via equações (2.64) e (2.65)"
            //"Step 10 - On the scale j, compute the coefficients of scale in scale j + 1 via equations (2.64) and (2.65)"
            nextScaleCoefficents[2 * (i)] = Math.pow(2, -(1 / 2)) * (lastScaleCoefficients[i] + waveletCoefficient);
            nextScaleCoefficents[2 * (i) + 1] = Math.pow(2, -(1 / 2)) * (lastScaleCoefficients[i] - waveletCoefficient);
        }

        return nextScaleCoefficents;
    }
    /*
    public static void main(String[] args) {

    double alpha = 0.9;
    double gamma = 0.13;
    double rho = 0.5;
    int n = 16;

    MMWGenerator mmw = new MMWGenerator(n, alpha, gamma, rho);

    Date inicio, fim, diferenca;
    SimpleDateFormat formato = new SimpleDateFormat("HH:mm:ss");

    long correcaoFuso = 10800000;

    for (int j = 0; j < 20; j++) {
    inicio = new Date();
    System.out.println(
    Resource.getString("BEGIN_GENERATION_AT")
    + formato.format(inicio));

    try {
    double[] traffic = mmw.getMMWTraffic();


    fim = new Date();
    System.out.println(
    Resource.getString("END_GENERATION_AT")
    + formato.format(fim));

    diferenca = new Date();
    diferenca.setTime(fim.getTime() - inicio.getTime() + correcaoFuso);
    System.out.println(
    Resource.getString("APPROXIMATED_EXECUTION_TIME")
    + formato.format(diferenca));
    System.out.println(
    Resource.getString("NUMBER_OF_SAMPLES_GENERATED")
    + (int) Math.pow(2, n));


    try {
    //TODO FIX: windows absolute path on code
    BufferedWriter out = new BufferedWriter(new FileWriter("C:/Documents and Settings/Tiago/Desktop/Tráfego/testeTrafego" + j + ".txt"));

    for (int i = 0; i < traffic.length; i++) {
    out.write(Double.toString(traffic[i]) + "\n");
    }

    out.close();

    } catch (IOException e) {
    }
    } catch (Exception e) {
    }
    }
    }*/
}
